function h = qqplot1(x,y,pvec)
%QQPLOT Display an empirical quantile-quantile plot.
%   modified on 04.03.2003 -cm
%   QQPLOT1(X) makes an empirical QQ-plot of the quantiles of
%   the data set X versus the quantiles of a standard Normal distribution.
%
%   QQPLOT1(X,Y) makes an empirical QQ-plot of the quantiles of
%   the data set X versus the quantiles of the data set Y.
%
%   H = QQPLOT1(X,Y,PVEC) allows you to specify the plotted quantiles in 
%   the vector PVEC. H is a handle to the plotted lines. 
%
%   When both X and Y are input, the default quantiles are those of the 
%   smaller data set.
%
%   The purpose of the quantile-quantile plot is to determine whether
%   the sample in X is drawn from a Normal (i.e., Gaussian) distribution,
%   or whether the samples in X and Y come from the same distribution
%   type.  If the samples do come from the same distribution (same shape),
%   even if one distribution is shifted and re-scaled from the other
%   (different location and scale parameters), the plot will be linear.

%   Copyright 1993-2002 The MathWorks, Inc. 
%   $Revision: 2.13 $  $Date: 2002/01/17 21:31:46 $

if nargin == 1
   y  =  sort(x);
   [x,n]  = plotpos(y);
   x  = norminv(x);
   xx = x;
   yy = y;
else
   n = -1;
   if nargin < 3
      nx = sum(~isnan(x));
      if (length(nx) > 1)
         nx = max(nx);
      end
      ny = sum(~isnan(y));
      if (length(ny) > 1)
         ny = max(ny);
      end
      n    = min(nx, ny);
      pvec = 100*((1:n) - 0.5) ./ n;
   end

   if (((size(x,1)==n) | (size(x,1)==1 & size(x,2)==n)) & ~any(isnan(x)))
      xx = sort(x);
   else
      xx = prctile(x,pvec);
   end
   if (((size(y,1)==n) | (size(y,1)==1 & size(y,2)==n)) & ~any(isnan(y)))
      yy = sort(y);
   else
      yy=prctile(y,pvec);
   end
end

q1x = prctile(x,25);
q3x = prctile(x,75);
q1y = prctile(y,25);
q3y = prctile(y,75);
qx = [q1x; q3x];
qy = [q1y; q3y];


dx = q3x - q1x;
dy = q3y - q1y;
slope = dy./dx;
centerx = (q1x + q3x)/2;
centery = (q1y + q3y)/2;
maxx = max(x);
minx = min(x);
maxy = centery + slope.*(maxx - centerx);
miny = centery - slope.*(centerx - minx);

mx = [minx; maxx];
my = [miny; maxy];


hh = plot(xx,yy,'+',qx,qy,'-',mx,my,'-.');
if nargout == 1
  h = hh;
end

if nargin == 1
   xlabel('Standard Normal Quantiles')
   ylabel('Quantiles of Input Sample')
   title ('QQ Plot of Sample Data versus Standard Normal')
else
   xlabel('X Quantiles');
   ylabel('Y Quantiles');
end

%===================== helper function 1====================
function [pp,n] = plotpos(sx)
%PLOTPOS Compute plotting positions for a probability plot
%   PP = PLOTPOS(SX) compute the plotting positions for a probabilty
%   plot of the columns of SX (or for SX itself if it is a vector).
%   SX must be sorted before being passed into PLOTPOS.  The ith
%   value of SX has plotting position (i-0.5)/n, where n is
%   the number of rows of SX.  NaN values are removed before
%   computing the plotting positions.
%
%   [PP,N] = PLOTPOS(SX) also returns N, the largest sample size
%   among the columns of SX.  This N can be used to set axis limits.

[n, m] = size(sx);
if n == 1
   sx = sx';
   n = m;
   m = 1;
end

nvec = sum(~isnan(sx));
pp = repmat((1:n)', 1, m);
pp = (pp-.5) ./ repmat(nvec, n, 1);
pp(isnan(sx)) = NaN;

if (nargout > 1)
   n = max(nvec);  % sample size for setting axis limits
end


%===================== helper function 2====================
function z = norminv(p,mu,sigma)
%NORMINV Inverse of the normal cumulative distribution function (cdf).
%   X = NORMINV(P,MU,SIGMA) finds the inverse of the normal cdf with
%   mean, MU, and standard deviation, SIGMA.
%
%   The size of X is the common size of the input arguments. A scalar input  
%   functions as a constant matrix of the same size as the other inputs.    
%
%   Default values for MU and SIGMA are 0 and 1 respectively.
%
%   See also NORMCDF, ERF, ERFC, ERFINV, ERFCINV.

%   References:
%      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
%      Functions", Government Printing Office, 1964, 7.1.1 and 26.2.2

%   Copyright 1993-2002 The MathWorks, Inc. 
%   $Revision: 2.12 $  $Date: 2002/01/17 21:31:31 $

if nargin < 2, mu = 0; end
if nargin < 3, sigma = 1; end

%[errorcode p mu sigma] = distchck(3,p,mu,sigma);
%if errorcode > 0
%    error('Requires non-scalar arguments to match in size.');
%end

% It is numerically preferable to use the complementary error function
% and norminv(p) = -sqrt(2)*erfcinv(2*p) to produce accurate results
% for p near zero.

z = (-sqrt(2)*sigma).*erfcinv(2*p) + mu;

%===================== helper function 3====================
function x = erfcinv(y)
%ERFCINV Inverse complementary error function.
%   X = ERFCINV(Y) is the inverse of the complementary error function
%   for each element of Y.  It satisfies y = erfc(x) 
%   for 2 >= y >= 0 and -Inf <= x <= Inf.
%
%   See also ERF, ERFC, ERFCX, ERFINV.

%   Copyright 1993-2002 The MathWorks, Inc. 
%   $Revision: 1.4 $  $Date: 2002/04/09 00:29:47 $

%   Original algorithm for norminv from Peter J. Acklam, jacklam@math.uio.no.

if ~isreal(y), error('Y must be real.'); end
x = zeros(size(y));

% Coefficients in rational approximations.
a = [  1.370600482778535e-02 -3.051415712357203e-01 ...
       1.524304069216834e+00 -3.057303267970988e+00  ...
       2.710410832036097e+00 -8.862269264526915e-01 ];
b = [ -5.319931523264068e-02  6.311946752267222e-01 ...
      -2.432796560310728e+00  4.175081992982483e+00 ...
      -3.320170388221430e+00 ];
c = [  5.504751339936943e-03  2.279687217114118e-01 ...
       1.697592457770869e+00  1.802933168781950e+00 ...
      -3.093354679843504e+00 -2.077595676404383e+00 ];
d = [  7.784695709041462e-03  3.224671290700398e-01 ...
       2.445134137142996e+00  3.754408661907416e+00 ];

% Define break-points.
ylow  = 0.0485;
yhigh = 1.9515;

% Rational approximation for central region
k = ylow <= y & y <= yhigh;
if any(k(:))
   q = y(k)-1;
   r = q.*q;
   x(k) = (((((a(1)*r+a(2)).*r+a(3)).*r+a(4)).*r+a(5)).*r+a(6)).*q ./ ...
          (((((b(1)*r+b(2)).*r+b(3)).*r+b(4)).*r+b(5)).*r+1);
end

% Rational approximation for lower region
k = 0 < y & y < ylow;
if any(k(:))
   q  = sqrt(-2*log(y(k)/2));
   x(k) = (((((c(1)*q+c(2)).*q+c(3)).*q+c(4)).*q+c(5)).*q+c(6)) ./ ...
           ((((d(1)*q+d(2)).*q+d(3)).*q+d(4)).*q+1);
end

% Rational approximation for upper region
k = yhigh < y & y < 2;
if any(k(:))
   q  = sqrt(-2*log(1-y(k)/2));
   x(k) = -(((((c(1)*q+c(2)).*q+c(3)).*q+c(4)).*q+c(5)).*q+c(6)) ./ ...
            ((((d(1)*q+d(2)).*q+d(3)).*q+d(4)).*q+1);
end

% The relative error of the approximation has absolute value less
% than 1.13e-9.  One iteration of Halley's rational method (third
% order) gives full machine precision.

% Newton's method: new x = x - f/f'
% Halley's method: new x = x - 1/(f'/f - (f"/f')/2)
% This function: f = erfc(x) - y, f' = -2/sqrt(pi)*exp(-x^2), f" = -2*x*f'

% Newton's correction
u = (erfc(x) - y) ./ (-2/sqrt(pi) * exp(-x.^2));

% Halley's step
x = x - u./(1+x.*u);

% Exceptional cases

x(y == 0) = Inf;
x(y == 2) = -Inf;
x(y < 0) = NaN;
x(y > 2) = NaN;
x(isnan(y)) = NaN;
